#ifndef VECTORMATH_H_
#define VECTORMATH_H_

#include <cmath>

namespace vm {

struct vec2;
struct vec3;
struct vec4;
struct mat2;
struct mat3;
struct mat4;

const float kPi = 3.1415926f;
const float kPiDiv2 = 0.5f * kPi;
const float kPiDiv180 = kPi / 180.0f;
const float k1DivPi = 1.0f / kPi;
const float k180DivPi = 180.0f / kPi;

//-----------------------------------------------------------------------------
// vec2
//-----------------------------------------------------------------------------
struct vec2 {
  vec2() {}
  explicit vec2(float s) : x(s), y(s) {}
  vec2(float x, float y) : x(x), y(y) {}
  explicit vec2(const vec3 &v);
  explicit vec2(const vec4 &v);

  const float &operator[](int i) const { return *(&x + i); }
  float &operator[](int i) { return *(&x + i); }

  vec2 operator-() const { return vec2(-x, -y); }

  vec2 &operator+=(float rhs) { x += rhs; y += rhs; return *this; }
  vec2 &operator+=(const vec2 &rhs) { x += rhs.x; y += rhs.y; return *this; }
  vec2 &operator-=(float rhs) { x -= rhs; y -= rhs; return *this; }
  vec2 &operator-=(const vec2 &rhs) { x -= rhs.x; y -= rhs.y; return *this; }
  vec2 &operator*=(float rhs) { x *= rhs; y *= rhs; return *this; }
  vec2 &operator*=(const vec2 &rhs) { x *= rhs.x; y *= rhs.y; return *this; }
  vec2 &operator/=(float rhs) {
    const float rcp = 1.0f / rhs;
    x *= rcp; y *= rcp; return *this;
  }
  vec2 &operator/=(const vec2 &rhs) { x /= rhs.x; y /= rhs.y; return *this; }

  float x, y;
};
inline bool operator==(const vec2 &lhs, const vec2 &rhs) {
  return lhs.x == rhs.x && lhs.y == rhs.y;
}
inline bool operator!=(const vec2 &lhs, const vec2 &rhs) {
  return !operator==(lhs, rhs);
}
inline vec2 operator+(float lhs, const vec2 &rhs) {
  return vec2(lhs + rhs.x, lhs + rhs.y);
}
inline vec2 operator+(const vec2 &lhs, float rhs) {
  return vec2(lhs.x + rhs, lhs.y + rhs);
}
inline vec2 operator+(const vec2 &lhs, const vec2 &rhs) {
  return vec2(lhs.x + rhs.x, lhs.y + rhs.y);
}
inline vec2 operator-(float lhs, const vec2 &rhs) {
  return vec2(lhs - rhs.x, lhs - rhs.y);
}
inline vec2 operator-(const vec2 &lhs, float rhs) {
  return vec2(lhs.x - rhs, lhs.y - rhs);
}
inline vec2 operator-(const vec2 &lhs, const vec2 &rhs) {
  return vec2(lhs.x - rhs.x, lhs.y - rhs.y);
}
inline vec2 operator*(float lhs, const vec2 &rhs) {
  return vec2(lhs * rhs.x, lhs * rhs.y);
}
inline vec2 operator*(const vec2 &lhs, float rhs) {
  return vec2(lhs.x * rhs, lhs.y * rhs);
}
inline vec2 operator*(const vec2 &lhs, const vec2 &rhs) {
  return vec2(lhs.x * rhs.x, lhs.y * rhs.y);
}
inline vec2 operator/(float lhs, const vec2 &rhs) {
  return vec2(lhs / rhs.x, lhs / rhs.y);
}
inline vec2 operator/(const vec2 &lhs, float rhs) {
  const float rcp = 1.0f / rhs;
  return vec2(lhs.x * rcp, lhs.y * rcp);
}
inline vec2 operator/(const vec2 &lhs, const vec2 &rhs) {
  return vec2(lhs.x / rhs.x, lhs.y / rhs.y);
}
//-----------------------------------------------------------------------------
// vec3
//-----------------------------------------------------------------------------
struct vec3 {
  vec3() {}
  explicit vec3(float s) : x(s), y(s), z(s) {}
  vec3(float x, float y, float z) : x(x), y(y), z(z) {}
  vec3(const vec2 &xy, float z) : x(xy[0]), y(xy[1]), z(z) {}
  vec3(float x, const vec2 &yz) : x(x), y(yz[0]), z(yz[1]) {}
  explicit vec3(const vec4 &v);

  const float &operator[](int i) const { return *(&x + i); }
  float &operator[](int i) { return *(&x + i); }

  vec3 operator-() const { return vec3(-x, -y, -z); }

  vec3 &operator+=(float rhs) {
    x += rhs; y += rhs; z += rhs; return *this;
  }
  vec3 &operator+=(const vec3 &rhs) {
    x += rhs.x; y += rhs.y; z += rhs.z; return *this;
  }
  vec3 &operator-=(float rhs) {
    x -= rhs; y -= rhs; z -= rhs; return *this;
  }
  vec3 &operator-=(const vec3 &rhs) {
    x -= rhs.x; y -= rhs.y; z -= rhs.z; return *this;
  }
  vec3 &operator*=(float rhs) {
    x *= rhs; y *= rhs; z *= rhs; return *this;
  }
  vec3 &operator*=(const vec3 &rhs) {
    x *= rhs.x; y *= rhs.y; z *= rhs.z; return *this;
  }
  vec3 &operator/=(float rhs) {
    const float rcp = 1.0f / rhs;
    x *= rcp; y *= rcp; z *= rcp; return *this;
  }
  vec3 &operator/=(const vec3 &rhs) {
    x /= rhs.x; y /= rhs.y; z /= rhs.z; return *this;
  }

  float x, y, z;
};
inline bool operator==(const vec3 &lhs, const vec3 &rhs) {
  return lhs.x == rhs.x && lhs.y == rhs.y && lhs.z == rhs.z;
}
inline bool operator!=(const vec3 &lhs, const vec3 &rhs) {
  return !operator==(lhs, rhs);
}
inline vec3 operator+(float lhs, const vec3 &rhs) {
  return vec3(lhs + rhs.x, lhs + rhs.y, lhs + rhs.z);
}
inline vec3 operator+(const vec3 &lhs, float rhs) {
  return vec3(lhs.x + rhs, lhs.y + rhs, lhs.z + rhs);
}
inline vec3 operator+(const vec3 &lhs, const vec3 &rhs) {
  return vec3(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z);
}
inline vec3 operator-(float lhs, const vec3 &rhs) {
  return vec3(lhs - rhs.x, lhs - rhs.y, lhs - rhs.z);
}
inline vec3 operator-(const vec3 &lhs, float rhs) {
  return vec3(lhs.x - rhs, lhs.y - rhs, lhs.z - rhs);
}
inline vec3 operator-(const vec3 &lhs, const vec3 &rhs) {
  return vec3(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z);
}
inline vec3 operator*(float lhs, const vec3 &rhs) {
  return vec3(lhs * rhs.x, lhs * rhs.y, lhs * rhs.z);
}
inline vec3 operator*(const vec3 &lhs, float rhs) {
  return vec3(lhs.x * rhs, lhs.y * rhs, lhs.z * rhs);
}
inline vec3 operator*(const vec3 &lhs, const vec3 &rhs) {
  return vec3(lhs.x * rhs.x, lhs.y * rhs.y, lhs.z * rhs.z);
}
inline vec3 operator/(float lhs, const vec3 &rhs) {
  return vec3(lhs / rhs.x, lhs / rhs.y, lhs / rhs.z);
}
inline vec3 operator/(const vec3 &lhs, float rhs) {
  const float rcp = 1.0f / rhs;
  return vec3(lhs.x * rcp, lhs.y * rcp, lhs.z * rcp);
}
inline vec3 operator/(const vec3 &lhs, const vec3 &rhs) {
  return vec3(lhs.x / rhs.x, lhs.y / rhs.y, lhs.z / rhs.z);
}
//-----------------------------------------------------------------------------
// vec4
//-----------------------------------------------------------------------------
struct vec4 {
  vec4() {}
  explicit vec4(float s) : x(s), y(s), z(s), w(s) {}
  vec4(float x, float y, float z, float w) : x(x), y(y), z(z), w(w) {}
  vec4(const vec2 &xy, float z, float w) : x(xy[0]), y(xy[1]), z(z), w(w) {}
  vec4(float x, const vec2 &yz, float w) : x(x), y(yz[0]), z(yz[1]), w(w) {}
  vec4(float x, float y, const vec2 &zw) : x(x), y(y), z(zw[0]), w(zw[1]) {}
  vec4(const vec3 &xyz, float w) : x(xyz[0]), y(xyz[1]), z(xyz[2]), w(w) {}
  vec4(float x, const vec3 &yzw) : x(x), y(yzw[0]), z(yzw[1]), w(yzw[2]) {}

  const float &operator[](int i) const { return *(&x + i); }
  float &operator[](int i) { return *(&x + i); }

  vec4 operator-() const { return vec4(-x, -y, -z, -w); }

  vec4 &operator+=(float rhs) {
    x += rhs; y += rhs; z += rhs; w += rhs; return *this;
  }
  vec4 &operator+=(const vec4 &rhs) {
    x += rhs.x; y += rhs.y; z += rhs.z; w += rhs.w; return *this;
  }
  vec4 &operator-=(float rhs) {
    x -= rhs; y -= rhs; z -= rhs; w -= rhs; return *this;
  }
  vec4 &operator-=(const vec4 &rhs) {
    x -= rhs.x; y -= rhs.y; z -= rhs.z; w -= rhs.w; return *this;
  }
  vec4 &operator*=(float rhs) {
    x *= rhs; y *= rhs; z *= rhs; w *= rhs; return *this;
  }
  vec4 &operator*=(const vec4 &rhs) {
    x *= rhs.x; y *= rhs.y; z *= rhs.z; w *= rhs.w; return *this;
  }
  vec4 &operator/=(float rhs) {
    const float rcp = 1.0f / rhs;
    x *= rcp; y *= rcp; z *= rcp; w *= rcp; return *this;
  }
  vec4 &operator/=(const vec4 &rhs) {
    x /= rhs.x; y /= rhs.y; z /= rhs.z; w /= rhs.w; return *this;
  }

  float x, y, z, w;
};
inline bool operator==(const vec4 &lhs, const vec4 &rhs) {
  return lhs.x == rhs.x && lhs.y == rhs.y && lhs.z == rhs.z && lhs.w == rhs.w;
}
inline bool operator!=(const vec4 &lhs, const vec4 &rhs) {
  return !operator==(lhs, rhs);
}
inline vec4 operator+(float lhs, const vec4 &rhs) {
  return vec4(lhs + rhs.x, lhs + rhs.y, lhs + rhs.z, lhs + rhs.w);
}
inline vec4 operator+(const vec4 &lhs, float rhs) {
  return vec4(lhs.x + rhs, lhs.y + rhs, lhs.z + rhs, lhs.w + rhs);
}
inline vec4 operator+(const vec4 &lhs, const vec4 &rhs) {
  return vec4(lhs.x + rhs.x, lhs.y + rhs.y, lhs.z + rhs.z, lhs.w + rhs.w);
}
inline vec4 operator-(float lhs, const vec4 &rhs) {
  return vec4(lhs - rhs.x, lhs - rhs.y, lhs - rhs.z, lhs - rhs.w);
}
inline vec4 operator-(const vec4 &lhs, float rhs) {
  return vec4(lhs.x - rhs, lhs.y - rhs, lhs.z - rhs, lhs.w - rhs);
}
inline vec4 operator-(const vec4 &lhs, const vec4 &rhs) {
  return vec4(lhs.x - rhs.x, lhs.y - rhs.y, lhs.z - rhs.z, lhs.w - rhs.w);
}
inline vec4 operator*(float lhs, const vec4 &rhs) {
  return vec4(lhs * rhs.x, lhs * rhs.y, lhs * rhs.z, lhs * rhs.w);
}
inline vec4 operator*(const vec4 &lhs, float rhs) {
  return vec4(lhs.x * rhs, lhs.y * rhs, lhs.z * rhs, lhs.w * rhs);
}
inline vec4 operator*(const vec4 &lhs, const vec4 &rhs) {
  return vec4(lhs.x * rhs.x, lhs.y * rhs.y, lhs.z * rhs.z, lhs.w * rhs.w);
}
inline vec4 operator/(float lhs, const vec4 &rhs) {
  return vec4(lhs / rhs.x, lhs / rhs.y, lhs / rhs.z, lhs / rhs.w);
}
inline vec4 operator/(const vec4 &lhs, float rhs) {
  const float rcp = 1.0f / rhs;
  return vec4(lhs.x * rcp, lhs.y * rcp, lhs.z * rcp, lhs.w * rcp);
}
inline vec4 operator/(const vec4 &lhs, const vec4 &rhs) {
  return vec4(lhs.x / rhs.x, lhs.y / rhs.y, lhs.z / rhs.z, lhs.w / rhs.w);
}

inline vec2::vec2(const vec3 &v) : x(v.x), y(v.y) {}
inline vec2::vec2(const vec4 &v) : x(v.x), y(v.y) {}
inline vec3::vec3(const vec4 &v) : x(v.x), y(v.y), z(v.z) {}
//------------------------------------------------------------------------------
// Angle and Trigonometry Functions
//------------------------------------------------------------------------------
inline float radians(float x) { return x * kPiDiv180; }
inline vec2 radians(const vec2 &x) {
  return vec2(radians(x[0]), radians(x[1]));
}
inline vec3 radians(const vec3 &x) {
  return vec3(radians(x[0]), radians(x[1]), radians(x[2]));
}
inline vec4 radians(const vec4 &x) {
  return vec4(radians(x[0]), radians(x[1]), radians(x[2]), radians(x[3]));
}

inline float degrees(float x) { return x * k180DivPi; }
inline vec2 degrees(const vec2 &x) {
  return vec2(degrees(x[0]), degrees(x[1]));
}
inline vec3 degrees(const vec3 &x) {
  return vec3(degrees(x[0]), degrees(x[1]), degrees(x[2]));
}
inline vec4 degrees(const vec4 &x) {
  return vec4(degrees(x[0]), degrees(x[1]), degrees(x[2]), degrees(x[3]));
}

inline float sin(float x) { return ::std::sin(x); }
inline vec2 sin(const vec2 &x) { return vec2(sin(x[0]), sin(x[1])); }
inline vec3 sin(const vec3 &x) { return vec3(sin(x[0]), sin(x[1]), sin(x[2])); }
inline vec4 sin(const vec4 &x) {
  return vec4(sin(x[0]), sin(x[1]), sin(x[2]), sin(x[3]));
}

inline float cos(float x) { return ::std::cos(x); }
inline vec2 cos(const vec2 &x) { return vec2(cos(x[0]), cos(x[1])); }
inline vec3 cos(const vec3 &x) { return vec3(cos(x[0]), cos(x[1]), cos(x[2])); }
inline vec4 cos(const vec4 &x) {
  return vec4(cos(x[0]), cos(x[1]), cos(x[2]), cos(x[3]));
}

inline float tan(float x) { return ::std::tan(x); }
inline vec2 tan(const vec2 &x) { return vec2(tan(x[0]), tan(x[1])); }
inline vec3 tan(const vec3 &x) { return vec3(tan(x[0]), tan(x[1]), tan(x[2])); }
inline vec4 tan(const vec4 &x) {
  return vec4(tan(x[0]), tan(x[1]), tan(x[2]), tan(x[3]));
}

inline float asin(float x) { return ::std::asin(x); }
inline vec2 asin(const vec2 &x) { return vec2(asin(x[0]), asin(x[1])); }
inline vec3 asin(const vec3 &x) {
  return vec3(asin(x[0]), asin(x[1]), asin(x[2]));
}
inline vec4 asin(const vec4 &x) {
  return vec4(asin(x[0]), asin(x[1]), asin(x[2]), asin(x[3]));
}

inline float acos(float x) { return ::std::acos(x); }
inline vec2 acos(const vec2 &x) { return vec2(acos(x[0]), acos(x[1])); }
inline vec3 acos(const vec3 &x) {
  return vec3(acos(x[0]), acos(x[1]), acos(x[2]));
}
inline vec4 acos(const vec4 &x) {
  return vec4(acos(x[0]), acos(x[1]), acos(x[2]), acos(x[3]));
}

inline float atan(float y, float x) { return ::std::atan2(y, x); }
inline vec2 atan(const vec2 &y, const vec2 &x) {
  return vec2(atan(y[0], x[0]), atan(y[1], x[1]));
}
inline vec3 atan(const vec2 &y, const vec3 &x) {
  return vec3(atan(y[0], x[0]), atan(y[1], x[1]), atan(y[2], x[2]));
}
inline vec4 atan(const vec2 &y, const vec4 &x) {
  return vec4(atan(y[0], x[0]), atan(y[1], x[1]), atan(y[2], x[2]),
              atan(y[3], x[3]));
}

inline float atan(float x) { return ::std::atan(x); }
inline vec2 atan(const vec2 &x) { return vec2(atan(x[0]), atan(x[1])); }
inline vec3 atan(const vec3 &x) {
  return vec3(atan(x[0]), atan(x[1]), atan(x[2]));
}
inline vec4 atan(const vec4 &x) {
  return vec4(atan(x[0]), atan(x[1]), atan(x[2]), atan(x[3]));
}
//-----------------------------------------------------------------------------
// Exponential Functions
//-----------------------------------------------------------------------------
inline float pow(float x, float y) { return ::std::pow(x, y); }
inline vec2 pow(const vec2 &x, const vec2 &y) {
  return vec2(pow(x[0], y[0]), pow(x[1], y[1]));
}
inline vec3 pow(const vec2 &x, const vec3 &y) {
  return vec3(pow(x[0], y[0]), pow(x[1], y[1]), pow(x[2], y[2]));
}
inline vec4 pow(const vec2 &x, const vec4 &y) {
  return vec4(pow(x[0], y[0]), pow(x[1], y[1]), pow(x[2], y[2]),
              pow(x[3], y[3]));
}

inline float exp(float x) { return ::std::exp(x); }
inline vec2 exp(const vec2 &x) { return vec2(exp(x[0]), exp(x[1])); }
inline vec3 exp(const vec3 &x) {
  return vec3(exp(x[0]), exp(x[1]), exp(x[2]));
}
inline vec4 exp(const vec4 &x) {
  return vec4(exp(x[0]), exp(x[1]), exp(x[2]), exp(x[3]));
}

inline float log(float x) { return ::std::log(x); }
inline vec2 log(const vec2 &x) { return vec2(log(x[0]), log(x[1])); }
inline vec3 log(const vec3 &x) {
  return vec3(log(x[0]), log(x[1]), log(x[2]));
}
inline vec4 log(const vec4 &x) {
  return vec4(log(x[0]), log(x[1]), log(x[2]), log(x[3]));
}

inline float exp2(float x) { return ::std::exp(0.693147180559945309f * x); }
inline vec2 exp2(const vec2 &x) { return vec2(exp2(x[0]), exp2(x[1])); }
inline vec3 exp2(const vec3 &x) {
  return vec3(exp2(x[0]), exp2(x[1]), exp2(x[2]));
}
inline vec4 exp2(const vec4 &x) {
  return vec4(exp2(x[0]), exp2(x[1]), exp2(x[2]), exp2(x[3]));
}

inline float log2(float x) { return ::std::log(x) / 0.693147180559945309f; }
inline vec2 log2(const vec2 &x) { return vec2(log2(x[0]), log2(x[1])); }
inline vec3 log2(const vec3 &x) {
  return vec3(log2(x[0]), log2(x[1]), log2(x[2]));
}
inline vec4 log2(const vec4 &x) {
  return vec4(log2(x[0]), log2(x[1]), log2(x[2]), log2(x[3]));
}

inline float sqrt(float x) { return ::std::sqrt(x); }
inline vec2 sqrt(const vec2 &x) { return vec2(sqrt(x[0]), sqrt(x[1])); }
inline vec3 sqrt(const vec3 &x) {
  return vec3(sqrt(x[0]), sqrt(x[1]), sqrt(x[2]));
}
inline vec4 sqrt(const vec4 &x) {
  return vec4(sqrt(x[0]), sqrt(x[1]), sqrt(x[2]), sqrt(x[3]));
}

inline float inversesqrt(float x) { return 1.0f / ::std::sqrt(x); }
inline vec2 inversesqrt(const vec2 &x) {
  return vec2(inversesqrt(x[0]), inversesqrt(x[1]));
}
inline vec3 inversesqrt(const vec3 &x) {
  return vec3(inversesqrt(x[0]), inversesqrt(x[1]), inversesqrt(x[2]));
}
inline vec4 inversesqrt(const vec4 &x) {
  return vec4(inversesqrt(x[0]), inversesqrt(x[1]), inversesqrt(x[2]),
              inversesqrt(x[3]));
}
//------------------------------------------------------------------------------
// Common Functions
//------------------------------------------------------------------------------
inline float abs(float x) { return x >= 0.0f ? x : -x; }
inline vec2 abs(const vec2 &x) { return vec2(abs(x[0]), abs(x[1])); }
inline vec3 abs(const vec3 &x) {
  return vec3(abs(x[0]), abs(x[1]), abs(x[2]));
}
inline vec4 abs(const vec4 &x) {
  return vec4(abs(x[0]), abs(x[1]), abs(x[2]), abs(x[3]));
}

inline float sign(float x) {
  if (x > 1.0f) return 1.0f; if (x < 1.0f) return -1.0f; return 0.0f;
}
inline vec2 sign(const vec2 &x) { return vec2(sign(x[0]), sign(x[1])); }
inline vec3 sign(const vec3 &x) {
  return vec3(sign(x[0]), sign(x[1]), sign(x[2]));
}
inline vec4 sign(const vec4 &x) {
  return vec4(sign(x[0]), sign(x[1]), sign(x[2]), sign(x[3]));
}

inline float floor(float x) { return ::std::floor(x); }
inline vec2 floor(const vec2 &x) { return vec2(floor(x[0]), floor(x[1])); }
inline vec3 floor(const vec3 &x) {
  return vec3(floor(x[0]), floor(x[1]), floor(x[2]));
}
inline vec4 floor(const vec4 &x) {
  return vec4(floor(x[0]), floor(x[1]), floor(x[2]), floor(x[3]));
}

inline float ceil(float x) { return ::std::ceil(x); }
inline vec2 ceil(const vec2 &x) { return vec2(ceil(x[0]), ceil(x[1])); }
inline vec3 ceil(const vec3 &x) {
  return vec3(ceil(x[0]), ceil(x[1]), ceil(x[2]));
}
inline vec4 ceil(const vec4 &x) {
  return vec4(ceil(x[0]), ceil(x[1]), ceil(x[2]), ceil(x[3]));
}

inline float fract(float x) { return x - floor(x); }
inline vec2 fract(const vec2 &x) { return vec2(fract(x[0]), fract(x[1])); }
inline vec3 fract(const vec3 &x) {
  return vec3(fract(x[0]), fract(x[1]), fract(x[2]));
}
inline vec4 fract(const vec4 &x) {
  return vec4(fract(x[0]), fract(x[1]), fract(x[2]), fract(x[3]));
}

inline float mod(float x, float y) {
  return x - y * floor(x / y);
}
inline vec2 mod(const vec2 &x, float y) {
  return vec2(mod(x[0], y), mod(x[1], y));
}
inline vec3 mod(const vec3 &x, float y) {
  return vec3(mod(x[0], y), mod(x[1], y), mod(x[2], y));
}
inline vec4 mod(const vec4 &x, float y) {
  return vec4(mod(x[0], y), mod(x[1], y), mod(x[2], y), mod(x[3], y));
}
inline vec2 mod(const vec2 &x, const vec2 &y) {
  return vec2(mod(x[0], y[0]), mod(x[1], y[1]));
}
inline vec3 mod(const vec3 &x, const vec3 &y) {
  return vec3(mod(x[0], y[0]), mod(x[1], y[1]), mod(x[2], y[2]));
}
inline vec4 mod(const vec4 &x, const vec4 &y) {
  return vec4(mod(x[0], y[0]), mod(x[1], y[1]), mod(x[2], y[2]),
              mod(x[3], y[3]));
}

inline float min(float x, float y) {
  return y < x ? y : x;
}
inline vec2 min(const vec2 &x, float y) {
  return vec2(min(x[0], y), min(x[1], y));
}
inline vec3 min(const vec3 &x, float y) {
  return vec3(min(x[0], y), min(x[1], y), min(x[2], y));
}
inline vec4 min(const vec4 &x, float y) {
  return vec4(min(x[0], y), min(x[1], y), min(x[2], y), min(x[3], y));
}
inline vec2 min(const vec2 &x, const vec2 &y) {
  return vec2(min(x[0], y[0]), min(x[1], y[1]));
}
inline vec3 min(const vec3 &x, const vec3 &y) {
  return vec3(min(x[0], y[0]), min(x[1], y[1]), min(x[2], y[2]));
}
inline vec4 min(const vec4 &x, const vec4 &y) {
  return vec4(min(x[0], y[0]), min(x[1], y[1]), min(x[2], y[2]),
              min(x[3], y[3]));
}

inline float max(float x, float y) {
  return x < y ? y : x;
}
inline vec2 max(const vec2 &x, float y) {
  return vec2(max(x[0], y), max(x[1], y));
}
inline vec3 max(const vec3 &x, float y) {
  return vec3(max(x[0], y), max(x[1], y), max(x[2], y));
}
inline vec4 max(const vec4 &x, float y) {
  return vec4(max(x[0], y), max(x[1], y), max(x[2], y), max(x[3], y));
}
inline vec2 max(const vec2 &x, const vec2 &y) {
  return vec2(max(x[0], y[0]), max(x[1], y[1]));
}
inline vec3 max(const vec3 &x, const vec3 &y) {
  return vec3(max(x[0], y[0]), max(x[1], y[1]), max(x[2], y[2]));
}
inline vec4 max(const vec4 &x, const vec4 &y) {
  return vec4(max(x[0], y[0]), max(x[1], y[1]), max(x[2], y[2]),
              max(x[3], y[3]));
}

inline float clamp(float x, float min_val, float max_val) {
  return min(max(x, min_val), max_val);
}
inline vec2 clamp(const vec2 &x, float min_val, float max_val) {
  return vec2(clamp(x[0], min_val, max_val), clamp(x[1], min_val, max_val));
}
inline vec3 clamp(const vec3 &x, float min_val, float max_val) {
  return vec3(clamp(x[0], min_val, max_val), clamp(x[1], min_val, max_val),
              clamp(x[2], min_val, max_val));
}
inline vec4 clamp(const vec4 &x, float min_val, float max_val) {
  return vec4(clamp(x[0], min_val, max_val), clamp(x[1], min_val, max_val),
              clamp(x[2], min_val, max_val), clamp(x[3], min_val, max_val));
}
inline vec2 clamp(const vec2 &x, const vec2 &min_val, const vec2 &max_val) {
  return vec2(clamp(x[0], min_val[0], max_val[0]),
              clamp(x[1], min_val[1], max_val[1]));
}
inline vec3 clamp(const vec3 &x, const vec3 &min_val, const vec3 &max_val) {
  return vec3(clamp(x[0], min_val[0], max_val[0]),
              clamp(x[1], min_val[1], max_val[1]),
              clamp(x[2], min_val[2], max_val[2]));
}
inline vec4 clamp(const vec4 &x, const vec4 &min_val, const vec4 &max_val) {
  return vec4(clamp(x[0], min_val[0], max_val[0]),
              clamp(x[1], min_val[1], max_val[1]),
              clamp(x[2], min_val[2], max_val[2]),
              clamp(x[3], min_val[3], max_val[3]));
}

inline float mix(float x, float y, float a) {
  return x * (1.0f - a) + y * a;
}
inline vec2 mix(const vec2 &x, const vec2 &y, float a) {
  return x * (1.0f - a) + y * a;
}
inline vec3 mix(const vec3 &x, const vec3 &y, float a) {
  return x * (1.0f - a) + y * a;
}
inline vec4 mix(const vec4 &x, const vec4 &y, float a) {
  return x * (1.0f - a) + y * a;
}
inline vec2 mix(const vec2 &x, const vec2 &y, const vec2 &a) {
  return vec2(mix(x[0], y[0], a[0]), mix(x[1], y[1], a[1]));
}
inline vec3 mix(const vec3 &x, const vec3 &y, const vec3 &a) {
  return vec3(mix(x[0], y[0], a[0]), mix(x[1], y[1], a[1]),
              mix(x[2], y[2], a[2]));
}
inline vec4 mix(const vec4 &x, const vec4 &y, const vec4 &a) {
  return vec4(mix(x[0], y[0], a[0]), mix(x[1], y[1], a[1]),
              mix(x[2], y[2], a[2]), mix(x[3], y[3], a[3]));
}

inline float step(float edge, float x) { return x < edge ? 0.0f : 1.0f; }
inline vec2 step(float edge, const vec2 &x) {
  return vec2(step(edge, x[0]), step(edge, x[1]));
}
inline vec3 step(float edge, const vec3 &x) {
  return vec3(step(edge, x[0]), step(edge, x[1]), step(edge, x[2]));
}
inline vec4 step(float edge, const vec4 &x) {
  return vec4(step(edge, x[0]), step(edge, x[1]), step(edge, x[2]),
              step(edge, x[3]));
}
inline vec2 step(const vec2 &edge, const vec2 &x) {
  return vec2(step(edge[0], x[0]), step(edge[1], x[1]));
}
inline vec3 step(const vec3 &edge, const vec3 &x) {
  return vec3(step(edge[0], x[0]), step(edge[1], x[1]), step(edge[2], x[2]));
}
inline vec4 step(const vec4 &edge, const vec4 &x) {
  return vec4(step(edge[0], x[0]), step(edge[1], x[1]), step(edge[2], x[2]),
              step(edge[3], x[3]));
}

inline float smoothstep(float edge0, float edge1, float x) {
  const float t = clamp((x - edge0) / (edge1 - edge0), 0.0f, 1.0f);
  return t * t * (3 - 2 * t);
}
inline vec2 smoothstep(float edge0, float edge1, const vec2 &x) {
  return vec2(smoothstep(edge0, edge1, x[0]), smoothstep(edge0, edge1, x[1]));
}
inline vec3 smoothstep(float edge0, float edge1, const vec3 &x) {
  return vec3(smoothstep(edge0, edge1, x[0]), smoothstep(edge0, edge1, x[1]),
              smoothstep(edge0, edge1, x[2]));
}
inline vec4 smoothstep(float edge0, float edge1, const vec4 &x) {
  return vec4(smoothstep(edge0, edge1, x[0]), smoothstep(edge0, edge1, x[1]),
              smoothstep(edge0, edge1, x[2]), smoothstep(edge0, edge1, x[3]));
}
inline vec2 smoothstep(const vec2 &edge0, const vec2 &edge1, const vec2 &x) {
  return vec2(smoothstep(edge0[0], edge1[0], x[0]),
              smoothstep(edge0[1], edge1[1], x[1]));
}
inline vec3 smoothstep(const vec3 &edge0, const vec3 &edge1, const vec3 &x) {
  return vec3(smoothstep(edge0[0], edge1[0], x[0]),
              smoothstep(edge0[1], edge1[1], x[1]),
              smoothstep(edge0[2], edge1[2], x[2]));
}
inline vec4 smoothstep(const vec4 &edge0, const vec4 &edge1, const vec4 &x) {
  return vec4(smoothstep(edge0[0], edge1[0], x[0]),
              smoothstep(edge0[1], edge1[1], x[1]),
              smoothstep(edge0[2], edge1[2], x[2]),
              smoothstep(edge0[3], edge1[3], x[3]));
}
//------------------------------------------------------------------------------
// Geometric Functions
//------------------------------------------------------------------------------
inline float length(float x) { return abs(x); }
inline float length(const vec2 &x) { return sqrt(x[0]*x[0] + x[1]*x[1]); }
inline float length(const vec3 &x) {
  return sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
}
inline float length(const vec4 &x) {
  return sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2] + x[3]*x[3]);
}

inline float distance(float x, float y) { return length(x - y); }
inline float distance(const vec2 &x, const vec2 &y) {
  return length(x - y);
}
inline float distance(const vec3 &x, const vec3 &y) {
  return length(x - y);
}
inline float distance(const vec4 &x, const vec4 &y) {
  return length(x - y);
}

inline float dot(float x, float y) { return x * y; }
inline float dot(const vec2 &x, const vec2 &y) {
  return x[0]*y[0] + x[1]*y[1];
}
inline float dot(const vec3 &x, const vec3 &y) {
  return x[0]*y[0] + x[1]*y[1] + x[2]*y[2];
}
inline float dot(const vec4 &x, const vec4 &y) {
  return x[0]*y[0] + x[1]*y[1] + x[2]*y[2] + x[3]*y[3];
}

inline vec3 cross(const vec3 &x, const vec3 &y) {
  return vec3(x[1]*y[2] - y[1]*x[2], x[2]*y[0] - y[2]*x[0],
              x[0]*y[1] - y[0]*x[1]);
}

inline float normalize(float x) { return x / length(x); }
inline vec2 normalize(const vec2 &x) { return x / length(x); }
inline vec3 normalize(const vec3 &x) { return x / length(x); }
inline vec4 normalize(const vec4 &x) { return x / length(x); }

inline float faceforward(float N, float l, float Nref) {
  return dot(Nref, l) < 0.0f ? N : -N;
}
inline vec2 faceforward(const vec2 &N, const vec2 &l, const vec2 &Nref) {
  return dot(Nref, l) < 0.0f ? N : -N;
}
inline vec3 faceforward(const vec3 &N, const vec3 &l, const vec3 &Nref) {
  return dot(Nref, l) < 0.0f ? N : -N;
}
inline vec4 faceforward(const vec4 &N, const vec4 &l, const vec4 &Nref) {
  return dot(Nref, l) < 0.0f ? N : -N;
}

inline float reflect(float l, float N) {
  return l - 2 * dot(N, l) * N;
}
inline vec2 reflect(const vec2 &l, const vec2 &N) {
  return l - 2 * dot(N, l) * N;
}
inline vec3 reflect(const vec3 &l, const vec3 &N) {
  return l - 2 * dot(N, l) * N;
}
inline vec4 reflect(const vec4 &l, const vec4 &N) {
  return l - 2 * dot(N, l) * N;
}

inline float refract(float l, float N, float eta) {
  const float k = 1.0f - eta*eta * (1.0f - dot(N, l) * dot(N, l));
  return k < 0.0f ? 0.0f : eta * l - (eta * dot(N, l) * sqrt(k)) * N;
}
inline vec2 refract(const vec2 &l, const vec2 &N, float eta) {
  const float k = 1.0f - eta*eta * (1.0f - dot(N, l) * dot(N, l));
  return k < 0.0f ? vec2(0.0f) : eta * l - (eta * dot(N, l) * sqrt(k)) * N;
}
inline vec3 refract(const vec3 &l, const vec3 &N, float eta) {
  const float k = 1.0f - eta*eta * (1.0f - dot(N, l) * dot(N, l));
  return k < 0.0f ? vec3(0.0f) : eta * l - (eta * dot(N, l) * sqrt(k)) * N;
}
inline vec4 refract(const vec4 &l, const vec4 &N, float eta) {
  const float k = 1.0f - eta*eta * (1.0f - dot(N, l) * dot(N, l));
  return k < 0.0f ? vec4(0.0f) : eta * l - (eta * dot(N, l) * sqrt(k)) * N;
}
//------------------------------------------------------------------------------
// mat2
//------------------------------------------------------------------------------
struct mat2 {
  mat2() {}
  explicit mat2(float s)
      : c0(s, 0.0f),
        c1(0.0f, s) {}
  mat2(const vec2 &c0, const vec2 &c1)
      : c0(c0), c1(c1) {}
  mat2(float c0x, float c0y,
       float c1x, float c1y)
      : c0(c0x, c0y),
        c1(c1x, c1y) {}

  const vec2 &operator[](int i) const { return *(&c0 + i); }
  vec2 &operator[](int i) { return *(&c0 + i); }

  mat2 operator-() const { return mat2(-c0, -c1); }

  mat2 &operator+=(float rhs) { c0 += rhs; c1 += rhs; return *this; }
  mat2 &operator+=(const mat2 &rhs) {
    c0 += rhs[0]; c1 += rhs[1]; return *this;
  }
  mat2 &operator-=(float rhs) { c0 -= rhs; c1 -= rhs; return *this; }
  mat2 &operator-=(const mat2 &rhs) {
    c0 -= rhs[0]; c1 -= rhs[1]; return *this;
  }
  mat2 &operator*=(float rhs) { c0 *= rhs; c1 *= rhs; return *this; }
  mat2 &operator*=(const mat2 &rhs) {
    mat2 operator*(const mat2 &m, const mat2 &n);
    return *this = *this * rhs;
  }
  mat2 &operator/=(float rhs) { c0 /= rhs; c1 /= rhs; return *this; }
  mat2 &operator/=(const mat2 &rhs) {
    c0 /= rhs[0]; c1 /= rhs[1]; return *this;
  }

  static const mat2 &Zero() { static mat2 m(0.0f); return m; }
  static const mat2 &Identity() { static mat2 m(1.0f); return m; }

  vec2 c0, c1;
};
inline bool operator==(const mat2 &lhs, const mat2 &rhs) {
  return lhs.c0 == rhs.c0 && lhs.c1 == rhs.c1;
}
inline bool operator!=(const mat2 &lhs, const mat2 &rhs) {
  return !operator==(lhs, rhs);
}
inline mat2 operator+(float lhs, const mat2 &rhs) {
  return mat2(lhs + rhs.c0, lhs + rhs.c1);
}
inline mat2 operator+(const mat2 &lhs, float rhs) {
  return mat2(lhs.c0 + rhs, lhs.c1 + rhs);
}
inline mat2 operator+(const mat2 &lhs, const mat2 &rhs) {
  return mat2(lhs.c0 + rhs.c0, lhs.c1 + rhs.c1);
}
inline mat2 operator-(float lhs, const mat2 &rhs) {
  return mat2(lhs - rhs.c0, lhs - rhs.c1);
}
inline mat2 operator-(const mat2 &lhs, float rhs) {
  return mat2(lhs.c0 - rhs, lhs.c1 - rhs);
}
inline mat2 operator-(const mat2 &lhs, const mat2 &rhs) {
  return mat2(lhs.c0 - rhs.c0, lhs.c1 - rhs.c1);
}
inline mat2 operator*(float lhs, const mat2 &rhs) {
  return mat2(lhs * rhs.c0, lhs * rhs.c1);
}
inline mat2 operator*(const mat2 &lhs, float rhs) {
  return mat2(lhs.c0 * rhs, lhs.c1 * rhs);
}
inline mat2 operator*(const mat2 &m, const mat2 &n) {
  mat2 r;
  r[0].x = m[0].x*n[0].x + m[1].x*n[0].y;
  r[1].x = m[0].x*n[1].x + m[1].x*n[1].y;
  r[0].y = m[0].y*n[0].x + m[1].y*n[0].y;
  r[1].y = m[0].y*n[1].x + m[1].y*n[1].y;
  return r;
}
inline vec2 operator*(const mat2 &m, const vec2 &v) {
  return vec2(m[0].x*v.x + m[1].x*v.y,
              m[0].y*v.x + m[1].y*v.y);
}
inline vec2 operator*(const vec2 &v, const mat2 &m) {
  return vec2(dot(v, m[0]), dot(v, m[1]));
}
inline mat2 operator/(float lhs, const mat2 &rhs) {
  return mat2(lhs / rhs.c0, lhs / rhs.c1);
}
inline mat2 operator/(const mat2 &lhs, float rhs) {
  return mat2(lhs.c0 / rhs, lhs.c1 / rhs);
}
inline mat2 operator/(const mat2 &lhs, const mat2 &rhs) {
  return mat2(lhs.c0 / rhs.c0, lhs.c1 / rhs.c1);
}
//------------------------------------------------------------------------------
// mat3
//------------------------------------------------------------------------------
struct mat3 {
  mat3() {}
  explicit mat3(float s)
      : c0(s, 0.0f, 0.0f),
        c1(0.0f, s, 0.0f),
        c2(0.0f, 0.0f, s) {}
  mat3(const vec3 &c0, const vec3 &c1, const vec3 &c2)
      : c0(c0), c1(c1), c2(c2) {}
  mat3(float c0x, float c0y, float c0z,
       float c1x, float c1y, float c1z,
       float c2x, float c2y, float c2z)
      : c0(c0x, c0y, c0z),
        c1(c1x, c1y, c1z),
        c2(c2x, c2y, c2z) {}

  const vec3 &operator[](int i) const { return *(&c0 + i); }
  vec3 &operator[](int i) { return *(&c0 + i); }

  mat3 operator-() const { return mat3(-c0, -c1, -c2); }

  mat3 &operator+=(float rhs) {
    c0 += rhs; c1 += rhs; c2 += rhs; return *this;
  }
  mat3 &operator+=(const mat3 &rhs) {
    c0 += rhs[0]; c1 += rhs[1]; c2 += rhs[2]; return *this;
  }
  mat3 &operator-=(float rhs) {
    c0 -= rhs; c1 -= rhs; c2 -= rhs; return *this;
  }
  mat3 &operator-=(const mat3 &rhs) {
    c0 -= rhs[0]; c1 -= rhs[1]; c2 -= rhs[2]; return *this;
  }
  mat3 &operator*=(float rhs) {
    c0 *= rhs; c1 *= rhs; c2 *= rhs; return *this;
  }
  mat3 &operator*=(const mat3 &rhs) {
    mat3 operator*(const mat3 &m, const mat3 &n);
    return *this = *this * rhs;
  }
  mat3 &operator/=(float rhs) {
    c0 /= rhs; c1 /= rhs; c2 /= rhs; return *this;
  }
  mat3 &operator/=(const mat3 &rhs) {
    c0 /= rhs[0]; c1 /= rhs[1]; c2 /= rhs[2]; return *this;
  }

  static const mat3 &Zero() { static mat3 m(0.0f); return m; }
  static const mat3 &Identity() { static mat3 m(1.0f); return m; }
  static mat3 Scale(float scax, float scay, float scaz) {
    return mat3(scax, 0.0f, 0.0f,
                0.0f, scay, 0.0f,
                0.0f, 0.0f, scaz);
  }
  static mat3 Scale(const vec3 &x) {
    return Scale(x[0], x[1], x[2]);
  }
  static mat3 RotationX(float x) {
    const float sinx = sin(x);
    const float cosx = cos(x);
    return mat3(1.0f,  0.0f, 0.0f,
                0.0f,  cosx, sinx,
                0.0f, -sinx, cosx);
  }
  static mat3 RotationY(float x) {
    const float sinx = sin(x);
    const float cosx = cos(x);
    return mat3(cosx, 0.0f, -sinx,
                0.0f, 1.0f,  0.0f,
                sinx, 0.0f,  cosx);
  }
  static mat3 RotationZ(float x) {
    const float sinx = sin(x);
    const float cosx = cos(x);
    return mat3( cosx, sinx, 0.0f,
                -sinx, cosx, 0.0f,
                 0.0f, 0.0f, 1.0f);
  }

  vec3 c0, c1, c2;
};
inline bool operator==(const mat3 &lhs, const mat3 &rhs) {
  return lhs.c0 == rhs.c0 && lhs.c1 == rhs.c1 && lhs.c2 == rhs.c2;
}
inline bool operator!=(const mat3 &lhs, const mat3 &rhs) {
  return !operator==(lhs, rhs);
}
inline mat3 operator+(float lhs, const mat3 &rhs) {
  return mat3(lhs + rhs.c0, lhs + rhs.c1, lhs + rhs.c2);
}
inline mat3 operator+(const mat3 &lhs, float rhs) {
  return mat3(lhs.c0 + rhs, lhs.c1 + rhs, lhs.c2 + rhs);
}
inline mat3 operator+(const mat3 &lhs, const mat3 &rhs) {
  return mat3(lhs.c0 + rhs.c0, lhs.c1 + rhs.c1, lhs.c2 + rhs.c2);
}
inline mat3 operator-(float lhs, const mat3 &rhs) {
  return mat3(lhs - rhs.c0, lhs - rhs.c1, lhs - rhs.c2);
}
inline mat3 operator-(const mat3 &lhs, float rhs) {
  return mat3(lhs.c0 - rhs, lhs.c1 - rhs, lhs.c2 - rhs);
}
inline mat3 operator-(const mat3 &lhs, const mat3 &rhs) {
  return mat3(lhs.c0 - rhs.c0, lhs.c1 - rhs.c1, lhs.c2 - rhs.c2);
}
inline mat3 operator*(float lhs, const mat3 &rhs) {
  return mat3(lhs * rhs.c0, lhs * rhs.c1, lhs * rhs.c2);
}
inline mat3 operator*(const mat3 &lhs, float rhs) {
  return mat3(lhs.c0 * rhs, lhs.c1 * rhs, lhs.c2 * rhs);
}
inline mat3 operator*(const mat3 &m, const mat3 &n) {
  mat3 r;
  r[0].x = m[0].x*n[0].x + m[1].x*n[0].y + m[2].x*n[0].z;
  r[1].x = m[0].x*n[1].x + m[1].x*n[1].y + m[2].x*n[1].z;
  r[2].x = m[0].x*n[2].x + m[1].x*n[2].y + m[2].x*n[2].z;
  r[0].y = m[0].y*n[0].x + m[1].y*n[0].y + m[2].y*n[0].z;
  r[1].y = m[0].y*n[1].x + m[1].y*n[1].y + m[2].y*n[1].z;
  r[2].y = m[0].y*n[2].x + m[1].y*n[2].y + m[2].y*n[2].z;
  r[0].z = m[0].z*n[0].x + m[1].z*n[0].y + m[2].z*n[0].z;
  r[1].z = m[0].z*n[1].x + m[1].z*n[1].y + m[2].z*n[1].z;
  r[2].z = m[0].z*n[2].x + m[1].z*n[2].y + m[2].z*n[2].z;
  return r;
}
inline vec3 operator*(const mat3 &m, const vec3 &v) {
  return vec3(m[0].x*v.x + m[1].x*v.y + m[2].x*v.z,
              m[0].y*v.x + m[1].y*v.y + m[2].y*v.z,
              m[0].z*v.x + m[1].z*v.y + m[2].z*v.z);
}
inline vec3 operator*(const vec3 &v, const mat3 &m) {
  return vec3(dot(v, m[0]), dot(v, m[1]), dot(v, m[2]));
}
inline mat3 operator/(float lhs, const mat3 &rhs) {
  return mat3(lhs / rhs.c0, lhs / rhs.c1, lhs / rhs.c2);
}
inline mat3 operator/(const mat3 &lhs, float rhs) {
  return mat3(lhs.c0 / rhs, lhs.c1 / rhs, lhs.c2 / rhs);
}
inline mat3 operator/(const mat3 &lhs, const mat3 &rhs) {
  return mat3(lhs.c0 / rhs.c0, lhs.c1 / rhs.c1, lhs.c2 / rhs.c2);
}
//------------------------------------------------------------------------------
// mat4
//------------------------------------------------------------------------------
struct mat4 {
  mat4() {}
  explicit mat4(float s)
      : c0(s, 0.0f, 0.0f, 0.0f),
        c1(0.0f, s, 0.0f, 0.0f),
        c2(0.0f, 0.0f, s, 0.0f),
        c3(0.0f, 0.0f, 0.0f, s) {}
  mat4(const vec4 &c0, const vec4 &c1, const vec4 &c2, const vec4 &c3)
      : c0(c0), c1(c1), c2(c2), c3(c3) {}
  mat4(float c0x, float c0y, float c0z, float c0w,
       float c1x, float c1y, float c1z, float c1w,
       float c2x, float c2y, float c2z, float c2w,
       float c3x, float c3y, float c3z, float c3w)
      : c0(c0x, c0y, c0z, c0w),
        c1(c1x, c1y, c1z, c1w),
        c2(c2x, c2y, c2z, c2w),
        c3(c3x, c3y, c3z, c3w) {}
  explicit mat4(const mat3 &m)
      : c0(m[0], 0.0f),
        c1(m[1], 0.0f),
        c2(m[2], 0.0f),
        c3(0.0f, 0.0f, 0.0f, 1.0f) {}

  const vec4 &operator[](int i) const { return *(&c0 + i); }
  vec4 &operator[](int i) { return *(&c0 + i); }

  mat4 operator-() const { return mat4(-c0, -c1, -c2, -c3); }

  mat4 &operator+=(float rhs) {
    c0 += rhs; c1 += rhs; c2 += rhs; c3 += rhs; return *this;
  }
  mat4 &operator+=(const mat4 &rhs) {
    c0 += rhs[0]; c1 += rhs[1]; c2 += rhs[2]; c3 += rhs[3]; return *this;
  }
  mat4 &operator-=(float rhs) {
    c0 -= rhs; c1 -= rhs; c2 -= rhs; c3 -= rhs; return *this;
  }
  mat4 &operator-=(const mat4 &rhs) {
    c0 -= rhs[0]; c1 -= rhs[1]; c2 -= rhs[2]; c3 -= rhs[3]; return *this;
  }
  mat4 &operator*=(float rhs) {
    c0 *= rhs; c1 *= rhs; c2 *= rhs; c3 *= rhs; return *this;
  }
  mat4 &operator*=(const mat4 &rhs) {
    mat4 operator*(const mat4 &m, const mat4 &n);
    return *this = *this * rhs;
  }
  mat4 &operator/=(float rhs) {
    c0 /= rhs; c1 /= rhs; c2 /= rhs; c3 /= rhs; return *this;
  }
  mat4 &operator/=(const mat4 &rhs) {
    c0 /= rhs[0]; c1 /= rhs[1]; c2 /= rhs[2]; c3 /= rhs[3]; return *this;
  }

  static const mat4 &Zero() { static mat4 m(0.0f); return m; }
  static const mat4 &Identity() { static mat4 m(1.0f); return m; }
  static mat4 Translation(float trax, float tray, float traz) {
    return mat4(1.0f, 0.0f, 0.0f, 0.0f,
                0.0f, 1.0f, 0.0f, 0.0f,
                0.0f, 0.0f, 1.0f, 0.0f,
                trax, tray, traz, 1.0f);
  }
  static mat4 Translation(const vec3 &x) {
    return Translation(x[0], x[1], x[2]);
  }
  static mat4 Scale(float x, float y, float z) {
    return mat4(mat3::Scale(x, y, z));
  }
  static mat4 Scale(const vec3 &x) { return Scale(x[0], x[1], x[2]); }
  static mat4 RotationX(float x) { return mat4(mat3::RotationX(x)); }
  static mat4 RotationY(float x) { return mat4(mat3::RotationY(x)); }
  static mat4 RotationZ(float x) { return mat4(mat3::RotationZ(x)); }
  static mat4 Perspective(float fovy, float aspect, float n, float f) {
    const float k = tan(kPiDiv2 - 0.5f * fovy);
    const float ri = 1.0f / (n - f);
    return mat4(k/aspect, 0.0f,        0.0f,  0.0f,
                    0.0f,    k,        0.0f,  0.0f,
                    0.0f, 0.0f,    (n+f)*ri, -1.0f,
                    0.0f, 0.0f, n*f*ri*2.0f,  0.0f);
  }
  static mat4 LookAt(const vec3 &eye, const vec3 &at, const vec3 &up) {
    const vec3 az = normalize(eye - at);
    const vec3 ax = normalize(cross(up, az));
    const vec3 ay = cross(az, ax);
    return mat4(        ax[0],         ay[0],         az[0], 0.0f,
                        ax[1],         ay[1],         az[1], 0.0f,
                        ax[2],         ay[2],         az[2], 0.0f,
                -dot(ax, eye), -dot(ay, eye), -dot(az, eye), 1.0f);
  }
  static mat4 LookAt(float ex, float ey, float ez,
                     float ax, float ay, float az,
                     float ux, float uy, float uz) {
    return LookAt(vec3(ex, ey, ez), vec3(ax, ay, az), vec3(ux, uy, uz));
  }

  vec4 c0, c1, c2, c3;
};
inline bool operator==(const mat4 &lhs, const mat4 &rhs) {
  return lhs.c0 == rhs.c0 && lhs.c1 == rhs.c1 && lhs.c2 == rhs.c2 &&
         lhs.c3 == rhs.c3;
}
inline bool operator!=(const mat4 &lhs, const mat4 &rhs) {
  return !operator==(lhs, rhs);
}
inline mat4 operator+(float lhs, const mat4 &rhs) {
  return mat4(lhs + rhs.c0, lhs + rhs.c1, lhs + rhs.c2, lhs + rhs.c3);
}
inline mat4 operator+(const mat4 &lhs, float rhs) {
  return mat4(lhs.c0 + rhs, lhs.c1 + rhs, lhs.c2 + rhs, lhs.c3 + rhs);
}
inline mat4 operator+(const mat4 &lhs, const mat4 &rhs) {
  return mat4(lhs.c0 + rhs.c0, lhs.c1 + rhs.c1, lhs.c2 + rhs.c2,
              lhs.c3 + rhs.c3);
}
inline mat4 operator-(float lhs, const mat4 &rhs) {
  return mat4(lhs - rhs.c0, lhs - rhs.c1, lhs - rhs.c2, lhs - rhs.c3);
}
inline mat4 operator-(const mat4 &lhs, float rhs) {
  return mat4(lhs.c0 - rhs, lhs.c1 - rhs, lhs.c2 - rhs, lhs.c3 - rhs);
}
inline mat4 operator-(const mat4 &lhs, const mat4 &rhs) {
  return mat4(lhs.c0 - rhs.c0, lhs.c1 - rhs.c1, lhs.c2 - rhs.c2,
              lhs.c3 - rhs.c3);
}
inline mat4 operator*(float lhs, const mat4 &rhs) {
  return mat4(lhs * rhs.c0, lhs * rhs.c1, lhs * rhs.c2, lhs * rhs.c3);
}
inline mat4 operator*(const mat4 &lhs, float rhs) {
  return mat4(lhs.c0 * rhs, lhs.c1 * rhs, lhs.c2 * rhs, lhs.c3 * rhs);
}
inline mat4 operator*(const mat4 &m, const mat4 &n) {
  mat4 r;
  r[0].x = m[0].x*n[0].x + m[1].x*n[0].y + m[2].x*n[0].z + m[3].x*n[0].w;
  r[1].x = m[0].x*n[1].x + m[1].x*n[1].y + m[2].x*n[1].z + m[3].x*n[1].w;
  r[2].x = m[0].x*n[2].x + m[1].x*n[2].y + m[2].x*n[2].z + m[3].x*n[2].w;
  r[3].x = m[0].x*n[3].x + m[1].x*n[3].y + m[2].x*n[3].z + m[3].x*n[3].w;
  r[0].y = m[0].y*n[0].x + m[1].y*n[0].y + m[2].y*n[0].z + m[3].y*n[0].w;
  r[1].y = m[0].y*n[1].x + m[1].y*n[1].y + m[2].y*n[1].z + m[3].y*n[1].w;
  r[2].y = m[0].y*n[2].x + m[1].y*n[2].y + m[2].y*n[2].z + m[3].y*n[2].w;
  r[3].y = m[0].y*n[3].x + m[1].y*n[3].y + m[2].y*n[3].z + m[3].y*n[3].w;
  r[0].z = m[0].z*n[0].x + m[1].z*n[0].y + m[2].z*n[0].z + m[3].z*n[0].w;
  r[1].z = m[0].z*n[1].x + m[1].z*n[1].y + m[2].z*n[1].z + m[3].z*n[1].w;
  r[2].z = m[0].z*n[2].x + m[1].z*n[2].y + m[2].z*n[2].z + m[3].z*n[2].w;
  r[3].z = m[0].z*n[3].x + m[1].z*n[3].y + m[2].z*n[3].z + m[3].z*n[3].w;
  r[0].w = m[0].w*n[0].x + m[1].w*n[0].y + m[2].w*n[0].z + m[3].w*n[0].w;
  r[1].w = m[0].w*n[1].x + m[1].w*n[1].y + m[2].w*n[1].z + m[3].w*n[1].w;
  r[2].w = m[0].w*n[2].x + m[1].w*n[2].y + m[2].w*n[2].z + m[3].w*n[2].w;
  r[3].w = m[0].w*n[3].x + m[1].w*n[3].y + m[2].w*n[3].z + m[3].w*n[3].w;
  return r;
}
inline vec4 operator*(const mat4 &m, const vec4 &v) {
  return vec4(m[0].x*v.x + m[1].x*v.y + m[2].x*v.z + m[3].x*v.w,
              m[0].y*v.x + m[1].y*v.y + m[2].y*v.z + m[3].y*v.w,
              m[0].z*v.x + m[1].z*v.y + m[2].z*v.z + m[3].z*v.w,
              m[0].w*v.x + m[1].w*v.y + m[2].w*v.z + m[3].w*v.w);
}
inline vec4 operator*(const vec4 &v, const mat4 &m) {
  return vec4(dot(v, m[0]), dot(v, m[1]), dot(v, m[2]), dot(v, m[3]));
}
inline mat4 operator/(float lhs, const mat4 &rhs) {
  return mat4(lhs / rhs.c0, lhs / rhs.c1, lhs / rhs.c2, lhs / rhs.c3);
}
inline mat4 operator/(const mat4 &lhs, float rhs) {
  return mat4(lhs.c0 / rhs, lhs.c1 / rhs, lhs.c2 / rhs, lhs.c3 / rhs);
}
inline mat4 operator/(const mat4 &lhs, const mat4 &rhs) {
  return mat4(lhs.c0 / rhs.c0, lhs.c1 / rhs.c1, lhs.c2 / rhs.c2,
              lhs.c3 / rhs.c3);
}
//------------------------------------------------------------------------------
// Matrix Functions
//------------------------------------------------------------------------------
inline mat2 matrixCompMult(const mat2 &x, const mat2 &y) {
  return mat2(x[0] * y[0], x[1] * y[1]);
}
inline mat3 matrixCompMult(const mat3 &x, const mat3 &y) {
  return mat3(x[0] * y[0], x[1] * y[1], x[2] * y[2]);
}
inline mat4 matrixCompMult(const mat4 &x, const mat4 &y) {
  return mat4(x[0] * y[0], x[1] * y[1], x[2] * y[2], x[3] * y[3]);
}
//------------------------------------------------------------------------------

} // namespace vm

#endif
